با استفاده از داده های زلزله ها در ایران و جهان به سوالات زیر پاسخ دهید.
توابع، داده ها و کتابخانه های زیر در طول تمرین مورد استفاده قرار می گیرند.
library(readr)
library(dplyr)
library(stringr)
library(ggmap)
library(highcharter)
library(plotly)
library(ggplot2)
library(rworldmap)
library(countrycode)
library(cowplot)
library(sp)
library(gganimate)
library(lubridate)
earthquakes = read_rds("data/historical_web_data_26112015.rds")
# convert coordinate to country
coords2country = function(points)
{
countriesSP = getMap(resolution='high')
pointsSP = SpatialPoints(points, proj4string=CRS(proj4string(countriesSP)))
indices = over(pointsSP, countriesSP)
indices$ISO3
}۱. با استفاده از داده های historical_web_data_26112015.rds و استفاده از نمودار پراکنش سه بعدی بسته plotly نمودار طول، عرض و عمق زلزله ها را رسم نمایید. علاوه بر آن بزرگی هر نقطه را برابر بزرگی زمین لرزه قرار دهید.
plot_ly(earthquakes,
x = ~Longitude, y = ~Latitude, z = ~Depth, size = ~Magnitude,
marker = list(color = ~Magnitude,
symbol = 'circle',
sizemode = 'diameter',
colorscale = c('#FFE8A8', '#884848'),
showscale = TRUE),
sizes = c(0.8, 8.8),
text = ~paste('Province: ', Province,
'<br>City: ', City,
'<br>Earthquake Magnitude: ', Magnitude)) %>%
add_markers() %>%
layout(scene = list(title = 'Iran Earthquakes',
xaxis = list(title = 'Longitude', range = c(40, 80)),
yaxis = list(title = 'Latitude', range = c(20, 40)),
zaxis = list(title = 'Depth', range = c(0, 180))))۲. پویانمایی سونامی های تاریخی را بر حسب شدت بر روی نقشه زمین رسم نمایید.(از داده زلزله های بزرگ استفاده نمایید.)
disaster = read_delim("data/disaster.txt",
"\t",
trim_ws = TRUE) %>%
select(lon = LONGITUDE,
lat = LATITUDE,
mag = INTENSITY,
dep = FOCAL_DEPTH,
country = COUNTRY,
year = YEAR,
t_death = TOTAL_DEATHS) %>%
na.omit()eq_an = ggplot() +
borders("world", colour="black", fill="gray") +
geom_point(aes(x = disaster$lon,
y = disaster$lat,
size = disaster$mag,
frame = disaster$mag),
color="blue") +
ggtitle("Earthquakes In The World") +
xlab("Longitute") +
ylab("Latitude")
gganimate(eq_an, "eq_an.gif") 
***
۳. نمودار چگالی دو بعدی زلزله های تاریخی ایران را رسم کنید.( از داده iran_earthquake.rds و لایه stat_density_2d استفاده نمایید).
eq_ir = read_rds("data/iran_earthquake.rds") %>%
filter(Long < 70)
ggplot(eq_ir, aes(x=Long, y=Lat)) +
geom_point() +
stat_density_2d(aes(fill = ..level.., colour = ..level..), geom = "polygon") +
ggtitle("IR Earthquakes Density") +
xlab("Longitude") +
ylab("Latitude")
۴. احتمال اینکه در ایران در پنج سال آینده زلزله به بزرگی هفت ریشتر رخ دهد را محاسبه کنید. (از احتمال شرطی استفاده کنید.)
برای این سوال احتمال اینکه در پنج سال آینده زلزله رخ دهد به شرط آن که بیشتر از ۷ ریشتر باشد را محاسبه می کنیم. برای اینکار احتمال اینکه در پنج سال آینده یک زلزله رخ دهد را به شرط آنکه آن زلزله بزرگ تر از ۷ ریشتر باشد را بایستی محاسبه کنیم. در پنج سال آینده باتوجه به داده ها با احتمال ۱۰۰٪ زلزله رخ خواهد داد پس بایستی فقط احتمال اینکه در پنج سال آینده زلزله ی به بزرگی بیشتر از ۷ ریشتر رخ دهد را محاسبه کنیم.
eq_ir_7 = disaster %>%
filter(mag >= 7, country == 'IRAN') %>%
arrange(-year) %>%
mutate(dis = lag(year),
dis = dis - year,
dis = ifelse(is.na(dis),2018-year,dis))
eq_ir_7_under5 = eq_ir_7 %>%
mutate(uder5 = ifelse(dis <= 5, 1, 0)) %>%
group_by(uder5) %>%
summarise(count = n())
# +3 for continues < 5 years
prob = eq_ir_7_under5 %>%
summarise(prob =((7 + 3)/sum(count)))
knitr::kable(prob, 'html')| prob |
|---|
| 0.625 |
۵. بر اساس داده های زلزله های بزرگ ابتدا تعداد و متوسط کشته زلزله ها را بر حسب کشور استخراج نمایید. سپس نمودار گرمایی تعداد کشته ها را بر روی کره زمین رسم نمایید.(مانند مثال زیر!)
disaster_country = disaster %>%
group_by(country) %>%
summarise(death_count = sum(as.integer(t_death)),death_mean = mean(as.integer(t_death))) %>%
arrange(-death_count) %>%
mutate(country_code = countrycode(country,
"country.name",
"iso3c"))
knitr::kable(disaster_country %>% select(-country_code))| country | death_count | death_mean |
|---|---|---|
| CHINA | 609347 | 18465.06061 |
| IRAN | 298903 | 15731.73684 |
| ARMENIA | 187000 | 62333.33333 |
| ITALY | 149711 | 7129.09524 |
| PAKISTAN | 146342 | 24390.33333 |
| TURKMENISTAN | 110001 | 55000.50000 |
| PERU | 83055 | 2373.00000 |
| AZERBAIJAN | 80087 | 26695.66667 |
| ECUADOR | 71202 | 17800.50000 |
| CHILE | 65469 | 2618.76000 |
| TURKEY | 59062 | 4218.71429 |
| INDIA | 55233 | 6904.12500 |
| VENEZUELA | 27808 | 3972.57143 |
| GUATEMALA | 23054 | 4610.80000 |
| TAJIKISTAN | 15529 | 5176.33333 |
| MOROCCO | 13728 | 6864.00000 |
| RUSSIA | 12025 | 1717.85714 |
| NEPAL | 11771 | 3923.66667 |
| PHILIPPINES | 10883 | 473.17391 |
| MEXICO | 10349 | 862.41667 |
| ALGERIA | 7308 | 1461.60000 |
| JAPAN | 6125 | 322.36842 |
| GUADELOUPE | 5000 | 5000.00000 |
| UZBEKISTAN | 4880 | 4880.00000 |
| GREECE | 3133 | 120.50000 |
| YEMEN | 2800 | 2800.00000 |
| TAIWAN | 2325 | 332.14286 |
| INDONESIA | 2202 | 100.09091 |
| COLOMBIA | 2158 | 269.75000 |
| EL SALVADOR | 1328 | 265.60000 |
| MACEDONIA | 1083 | 361.00000 |
| AFGHANISTAN | 934 | 155.66667 |
| USA | 655 | 27.29167 |
| ARGENTINA | 516 | 103.20000 |
| KAZAKHSTAN | 451 | 225.50000 |
| GUINEA | 443 | 443.00000 |
| GEORGIA | 280 | 93.33333 |
| NEW ZEALAND | 258 | 86.00000 |
| ALBANIA | 243 | 121.50000 |
| PAPUA NEW GUINEA | 180 | 30.00000 |
| USA TERRITORY | 168 | 84.00000 |
| MONTENEGRO | 131 | 131.00000 |
| BULGARIA | 130 | 65.00000 |
| PORTUGAL | 113 | 113.00000 |
| COSTA RICA | 93 | 31.00000 |
| KYRGYZSTAN | 79 | 39.50000 |
| ETHIOPIA | 70 | 35.00000 |
| AZORES (PORTUGAL) | 69 | 69.00000 |
| SOLOMON ISLANDS | 62 | 20.66667 |
| MONGOLIA | 30 | 30.00000 |
| CANADA | 28 | 28.00000 |
| GHANA | 22 | 22.00000 |
| SLOVENIA | 22 | 11.00000 |
| IRAQ | 20 | 20.00000 |
| ROMANIA | 18 | 6.00000 |
| AUSTRALIA | 12 | 12.00000 |
| EGYPT | 12 | 12.00000 |
| UKRAINE | 11 | 11.00000 |
| DOMINICAN REPUBLIC | 8 | 4.00000 |
| CROATIA | 5 | 2.50000 |
| HONDURAS | 3 | 3.00000 |
| MYANMAR (BURMA) | 3 | 3.00000 |
| SERBIA | 3 | 1.50000 |
| BELGIUM | 2 | 2.00000 |
| CYPRUS | 2 | 2.00000 |
| DJIBOUTI | 2 | 2.00000 |
| AUSTRIA | 1 | 1.00000 |
| BOSNIA-HERZEGOVINA | 1 | 1.00000 |
| NETHERLANDS | 1 | 1.00000 |
| SOUTH AFRICA | 1 | 1.00000 |
| THAILAND | 1 | 1.00000 |
mapCountryData(joinCountryData2Map(disaster_country, nameJoinColumn = "country_code"),
nameColumnToPlot = "death_count",
mapTitle="Earthquakes Total Death")70 codes from your data successfully matched countries in the map
1 codes from your data failed to match with a country code in the map
175 codes from the map weren't represented in your data
mapCountryData(joinCountryData2Map(disaster_country, nameJoinColumn = "country_code"),
nameColumnToPlot = "death_mean",
mapTitle="Earthquakes Death Rate")70 codes from your data successfully matched countries in the map
1 codes from your data failed to match with a country code in the map
175 codes from the map weren't represented in your data
۶. با استفاده از داده لرزه های بزرگ و به وسیله طول، عرض، شدت، عمق مدلی برای پیش بینی تعداد کشته های زلزله بیابید.
مدل glm را باتوجه به نوع داده ها با خانواده ی گاما ایجاد می کنیم و آن را ارزیابی می کنیم و درمی یابیم که به عمق و شدت آن بستگی دارد و مکان آن اهمیتی در مدل ندارد پس مدل جدید را براین مبنا مدل سازی می کنیم و ارزیابی می کنیم که نتیجه ی قابل قبولی به دست می دهد.
eq_glm = glm(t_death ~ mag + dep + lat + lon, family = Gamma(link = "inverse"), data = disaster)
summary(eq_glm)
Call:
glm(formula = t_death ~ mag + dep + lat + lon, family = Gamma(link = "inverse"),
data = disaster)
Deviance Residuals:
Min 1Q Median 3Q Max
-4.055 -3.273 -2.611 -1.587 7.681
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.158e-03 2.561e-04 4.524 7.90e-06 ***
mag -1.070e-04 2.406e-05 -4.445 1.12e-05 ***
dep 7.511e-06 3.933e-06 1.910 0.0568 .
lat -9.875e-07 2.824e-06 -0.350 0.7267
lon -2.387e-07 4.959e-07 -0.481 0.6304
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for Gamma family taken to be 14.9034)
Null deviance: 3965.6 on 429 degrees of freedom
Residual deviance: 3398.8 on 425 degrees of freedom
AIC: 5834
Number of Fisher Scoring iterations: 10eq_glm = glm(t_death ~ mag + dep, family = Gamma(link = "inverse"), data = disaster)
summary(eq_glm)
Call:
glm(formula = t_death ~ mag + dep, family = Gamma(link = "inverse"),
data = disaster)
Deviance Residuals:
Min 1Q Median 3Q Max
-3.895 -3.270 -2.619 -1.607 7.898
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.223e-03 2.326e-04 5.259 2.30e-07 ***
mag -1.153e-04 2.197e-05 -5.247 2.44e-07 ***
dep 6.468e-06 1.967e-06 3.287 0.0011 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for Gamma family taken to be 15.12215)
Null deviance: 3965.6 on 429 degrees of freedom
Residual deviance: 3409.4 on 427 degrees of freedom
AIC: 5832.1
Number of Fisher Scoring iterations: 10۷. با استفاده از داده worldwide.csv به چند سوال زیر پاسخ دهید. تحقیق کنید آیا می توان از پیش لرزه، زلزله اصلی را پیش بینی کرد؟
از روی داده ها پیش لرزه ها را استخراج کرده و روی آن ها تعداد و میانگین آن ها و قدرت زلزله ی اصلی را درکنار هم جمع کرده و مدلی را استخراج می کنیم که این پیش بینی را صورت دهد که خروجی آن نتیجه ی قابل قبولی را گزارش می کند.
worldwide = read_csv("data/worldwide.csv")
worldwide = worldwide %>%
mutate(year = year(time),
month = month(time),
day = day(time),
location = as.data.frame(str_split_fixed(place, ',',2))$V2,
location = ifelse(location == '' , as.character(place), as.character(location))) %>%
select(time, year, month, day, latitude, longitude, depth, mag, location) %>%
group_by(year, month, location) %>%
arrange(time) %>%
mutate(mag_max = max(mag))
worldwide$groupid = worldwide %>% group_indices(year, month, location)
worldwide_temp = worldwide %>%
filter(mag == mag_max) %>%
ungroup() %>%
select(groupid,time_max = time)
worldwide = merge(worldwide, worldwide_temp)
worldwide = worldwide %>%
mutate(pre_past = ifelse(time > time_max, 'pre',ifelse(time == time_max, 'eq', 'past')))
worldwide_pred = worldwide %>%
filter(pre_past == 'pre' | pre_past == 'eq') %>%
ungroup() %>%
select(groupid, mag, pre_past) %>%
group_by(groupid,pre_past) %>%
summarise(avg_mag = mean(mag), count = n()) %>%
ungroup()
wwpm = worldwide_pred %>%
filter(pre_past == 'pre') %>%
select(groupid, count, avg_mag)
wwpmm = worldwide_pred %>%
filter(pre_past == 'eq') %>%
select(groupid, main = avg_mag)
x = merge(wwpm,wwpmm) %>%
select(-groupid)
glm_eq_pred = glm(count~., data = x)
summary(glm_eq_pred)
Call:
glm(formula = count ~ ., data = x)
Deviance Residuals:
Min 1Q Median 3Q Max
-92.62 -15.38 -1.20 10.15 1100.56
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 14.463 6.533 2.214 0.0269 *
avg_mag -44.653 2.148 -20.790 <2e-16 ***
main 37.713 1.681 22.429 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for gaussian family taken to be 2502.604)
Null deviance: 6924258 on 2218 degrees of freedom
Residual deviance: 5545771 on 2216 degrees of freedom
AIC: 23666
Number of Fisher Scoring iterations: 2۸. گزاره " آیا شدت زلزله به عمق آن بستگی دارد" را تحقیق کنید؟ (طبیعتا از آزمون فرض باید استفاده کنید.)
cor.test(worldwide$depth, worldwide$mag, method = 'spearman')
Spearman's rank correlation rho
data: worldwide$depth and worldwide$mag
S = 4.176e+13, p-value < 2.2e-16
alternative hypothesis: true rho is not equal to 0
sample estimates:
rho
0.2376803 mat_eq = as.matrix(worldwide %>%
select(depth, mag) %>%
mutate(depth = depth + abs(min(depth)))
)
chisq.test(mat_eq)
Pearson's Chi-squared test
data: mat_eq
X-squared = 489570, df = 69011, p-value < 2.2e-16تست وابستگی بیان می کند که غیرمرتبط بودن این دو متغیر از یکدیگر رد می شود پس می توان نتیجه گرفت که بین آن ها ارتباط وجود دارد و تست استقلال نیز بیان می کند که این دو متغیر نیز مستقل بودنشان از هم رد می شود و وابستگی آن ها نیز تایید می شود. پس فرض مستقل بودن این دو متغیر از هم رد می شود.
۹. میانگین سالانه زلزله ها را بر حسب کشور به دست آورید. آیا میتوان دلیلی در تایید یا رد تئوری هارپ ارائه کرد.
long_lat = worldwide %>%
select(longitude, latitude)
long_lat$country = coords2country(long_lat)
worldwide$country_code = long_lat$country
country_code = worldwide %>%
ungroup() %>%
group_by(location) %>%
select(location, country_code) %>%
na.omit() %>%
distinct(country_code, location) %>%
ungroup() %>%
select(location, c_code = country_code) %>%
group_by(location) %>%
slice(1) %>%
ungroup()
worldwide_fixed_country_code = merge(worldwide, country_code)
worldwide_fixed_country_code = worldwide_fixed_country_code %>%
mutate(c_code = ifelse(!is.na(country_code),
as.character(country_code),
as.character(c_code)
),
country_code = c_code
) %>%
select(-c_code)
worldwide = worldwide_fixed_country_code
worldwide_accu = worldwide %>%
group_by(year, country_code) %>%
summarise(eq_count = n(),
mean_mag = mean(mag),
mean_dep = mean(depth))
eq_2015 = worldwide_accu %>%
filter(year == 2015) %>%
select(eq_count)
eq_2016 = worldwide_accu %>%
filter(year == 2016) %>%
select(eq_count)
eq_2017 = worldwide_accu %>%
filter(year == 2017) %>%
select(eq_count)
eq_2018 = worldwide_accu %>%
filter(year == 2018) %>%
select(eq_count)
t.test(eq_2015$eq_count,eq_2016$eq_count)
Welch Two Sample t-test
data: eq_2015$eq_count and eq_2016$eq_count
t = -2.462, df = 109.01, p-value = 0.01538
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
-381.34944 -41.19174
sample estimates:
mean of x mean of y
32.98413 244.25472 t.test(eq_2015$eq_count,eq_2017$eq_count)
Welch Two Sample t-test
data: eq_2015$eq_count and eq_2017$eq_count
t = -2.1236, df = 108.3, p-value = 0.03598
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
-387.70921 -13.36027
sample estimates:
mean of x mean of y
32.98413 233.51887 t.test(eq_2015$eq_count,eq_2018$eq_count)
Welch Two Sample t-test
data: eq_2015$eq_count and eq_2018$eq_count
t = -1.2653, df = 96.949, p-value = 0.2088
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
-214.29194 47.43411
sample estimates:
mean of x mean of y
32.98413 116.41304 با توجه به اینکه پروژه ی هارپ در سال ۲۰۱۵ متوقف شد زلزله های این سال را با سال های بعدی آن مقایسه کردم و تست t را بر روی آن ها اجرا کردم تا تفاوت این سال ها را باهم مقایسه کنم. نتایج این تست نشان داد که این سال ها بایکدیگر تفاوت معناداری ندارند و این تئوری را می توان رد کرد.
۱۰. سه حقیقت جالب در مورد زلزله بیابید.
۱۰.۱. نمودار جعبه ای زلزله ها در سال ها و جایگاه ایران در آن
ww_mag_year = worldwide %>%
group_by(year, country_code) %>%
summarise(mean_mag = mean(mag)) %>%
select(year, mean_mag) %>%
ungroup()
iran_w = worldwide %>%
group_by(year, country_code) %>%
summarise(mean_mag = mean(mag)) %>%
filter(country_code == 'IRN') %>%
select(year, mean_mag) %>%
ungroup()
ggplot(ww_mag_year,aes(x = year, y = mean_mag, group = year)) +
geom_boxplot() +
geom_line(data = iran_w,aes(x = year, y = mean_mag, group = 1, color = 'Iran'), size = 1.5) +
ggtitle("Average Earthquake Magnitude In Worldwide") +
ylab("Magnitude Mean") +
xlab("Year")
۱۰.۲. ده زلزله خیزترین کشورها
country_reg = country_code %>%
group_by(c_code) %>%
summarise(reg_num = n()) %>%
ungroup() %>%
select(country_code = c_code, reg_num)
ww_mag_count = worldwide %>%
group_by(country_code) %>%
summarise(mean_mag = mean(mag), count = n())
ww_mag_count = merge(ww_mag_count, country_reg)
top_10 = ww_mag_count %>%
mutate(score = (count*mean_mag)/reg_num) %>%
arrange(-score) %>%
slice(1:10)
knitr::kable(top_10, "html")| country_code | mean_mag | count | reg_num | score |
|---|---|---|---|---|
| JPN | 4.531295 | 2486 | 1 | 11264.80 |
| PNG | 4.576130 | 2191 | 1 | 10026.30 |
| NZL | 4.537366 | 2058 | 1 | 9337.90 |
| PRI | 2.886450 | 2972 | 1 | 8578.53 |
| PHL | 4.531641 | 1438 | 1 | 6516.50 |
| IDN | 4.462053 | 4277 | 3 | 6361.40 |
| TON | 4.628141 | 1329 | 1 | 6150.80 |
| VGB | 2.986832 | 1632 | 1 | 4874.51 |
| VUT | 4.580264 | 1059 | 1 | 4850.50 |
| RUS | 4.425916 | 1092 | 1 | 4833.10 |
۱۰.۳. نحسی عدد ۱۳ در وقوع زلزله در روز ۱۳ام.
thirteen = worldwide %>%
filter(day == 13)
not_thirteen = worldwide %>%
filter(day != 13)
t.test(thirteen$mag, not_thirteen$mag, alternative = "greater")
Welch Two Sample t-test
data: thirteen$mag and not_thirteen$mag
t = 6.5195, df = 2091, p-value = 4.407e-11
alternative hypothesis: true difference in means is greater than 0
95 percent confidence interval:
0.0961231 Inf
sample estimates:
mean of x mean of y
3.873674 3.745097 این تست رد می شود و نمی توان نتیجه گرفت که روز ۱۳ام روز نحسی برای زلزله است.